
N67- 1 1 7 1 


■ ..;,'t^ : v; 
'?'£&¥;■ 'rL... ■ " r '-.;.:^‘ 

- k „• ^ .>- ■ i - •• ;v i - 


NASA ' 
Technical 


261 3 


December T 986 


Method for Inviscid 


Suresh M; Deshpande 


3rv* -■ 
/i i- < V • 


{NASfi-lf-2613) A SECCKB-CiCEE ACCUPAGE 
K JKETIC-lHECtil-i.A£E£ fETECE JCf JAVI3CIC 
CCMPH E S3 1 ELf ilCh£ (NASA) 42 f C3CI 20D 


HI/ 34 


Unclat 
4 3 7 Id 


NASA 

Technical 

Paper 

2613 

1986 


NASA 

National Aeronautics 
and Space Administration 


A Second-Order Accurate 
Kinetic-Theory-Based 
Method for Inviscid 
Compressible Flows 


Suresh M. Deshpande 

Langley Research Center 
Hampton, Virginia 


Scientific and Technical 
Information Branch 


Summary 

A new upwind method for the numerical solution 
of the Euler equations is presented. This method, 
called the kinetic numerical method (KNM), is based 
on the well-known fact that the Euler equations are 
the moments of the Boltzmann equation of the ki- 
netic theory of gases when the distribution function 
is Maxwellian. The KNM consists of two phases, the 
convection phase and the collision phase. The veloc- 
ity distribution function at the end of the convection 
phase is the solution of the collisionless Boltzmann 
equation, which is linear and hyperbolic. The col- 
lision phase instantaneously relaxes the distribution 
to the local Maxwellian distribution. The fluid dy- 
namic variables of density, velocity, and internal en- 
ergy are obtained as moments of the velocity distri- 
bution function at the end of the convection phase. 
The KNM is an explicit method and is uncondition- 
ally stable. It is an upwind method because of the 
way the solution is obtained in the convection phase 
and satisfies the entropy condition due to the con- 
vexity of the Boltzmann H- function, which decreases 
as the velocity distribution function suddenly relaxes 
to the local Maxwellian distribution during the colli- 
sion phase. Because of its explicit nature, the KNM 
is highly vectorizable and can be easily made into 
a total variation diminishing (TVD) method for the 
distribution function through a suitable choice of the 
interpolation strategy. 

A fascinating aspect of the present work is the 
use of the antidiffusive Chapman-Enskog theory in 
the development of the second-order accurate KNM. 
Normally the Chapman-Enskog theory is associated 
with the Navier-Stokes equations, and use of the 
Chapman-Enskog theory for constructing the second- 
order KNM is very unusual. It is shown that to cancel 
the large amount of viscosity in the first-order KNM, 
antidiffusive terms are required, and these can be in- 
troduced through the Chapman-Enskog distribution 
if the shear stress tensor r and heat flux vector q are 
made antidiffusive. Through the addition of diffusive 
terms corresponding to viscosity and heat conduction 
to the expressions for r and q used in the Chapman- 
Enskog distribution, the KNM can be easily extended 
to the solution of the Navier-Stokes equations. 

The KNM is applied to a one-dimensional shock- 
propagation problem, and the results demonstrate 
the capability of the method for giving “wiggle-free” 
accurate solutions. Also, the Euler space-marching 
KNM is applied to obtain the steady solution to a 
shock-reflection problem. In this application, the 
boundary condition is treated differently from other 
finite-difference methods in that the flow tangency 
condition at the flat plate is imposed by reversing 


the normal velocity of the reflected molecules. The 
boundary condition is at the level of the velocity 
distribution function, which is natural in the KNM. 
The pressure profiles are compared with a second- 
order and a third-order accurate TVD method, and 
the ability of the KNM to yield an accurate wiggle- 
free solution is demonstrated. 


Introduction 

Recently, several numerical methods based on ki- 
netic theory have been proposed for the computation 
of inviscid compressible flows (refs. 1 to 4). These 
methods use the Boltzmann or the Boltzmann-like 
equation of the kinetic theory of gases. These ki- 
netic numerical methods (KNM’s), as Aristov and 
Tcheremissine (ref. 4) call them, are based on the 
connection between the Boltzmann equation and the 
Euler equations governing the dynamics of inviscid 
compressible flows. The basic unknown in the Boltz- 
mann equation is the velocity distribution function 
/(t,x, v,/), where t is time, x is the position vec- 
tor, v is the molecular velocity vector, and I is the 
independent internal-energy variable corresponding 
to nontranslational degrees of freedom. The spatio- 
temporal evolution of / is governed by the Boltzmann 
equation 

Equation (1) consists of the streaming, or convective, 
term + v • / and the collision term J(/, /). 

The convective term gives the rate of change of / per 
unit volume in (x, v,/) space because of movement 
of molecules, and the collision term gives the rate of 
change of / because of intermolecular collisions. 

The field variables p, u, and e (which are, re- 
spectively, the mass density, the fluid velocity vector, 
and the specific total internal energy) are related to 
/ through moment equations (see appendix A). It is 
shown in appendix A that the moments of the Boltz- 
mann equation (1) are the Euler equations if the ve- 
locity distribution function / is the Maxwellian dis- 
tribution F as follows: 


(y ~ u) 2 exp (— // 1 0 ) 

2 RT j I 0 

( 2 ) 

where R is the gas constant, T is the temperature, 
I is the internal-energy variable corresponding to 
nontranslational degrees of freedom, 


F = 


(2? xRTf/ 2 


exp 


Io 


(2 + Df) ~ iDf 
2( 7 -l) 


RT 



D j is the degrees of freedom of a molecule, and 7 is 
the ratio of specific heats. It is interesting to note 
here that if / is assumed to be the Chapman-Enskog 
distribution (ref. 5), then the moment equations re- 
duce to the Navier-Stokes equations. This connection 
between the Boltzmann equation and the equations 
of fluid dynamics is the basis of all KNM’s. 

Several variants of KNM’s can be obtained by 
adopting particle-in-cell (PIC), fluid-in-cell (FLIC), 
or finite-difference methods. Pullin (ref. 1) follows 
the statistical PIC method, in which the computa- 
tional domain is divided into several Eulerian cells. 
The fluid is represented in the PIC method by a 
large number of discrete Lagrangian mass points, 
called particles, which move through the Eulerian 
mesh convecting mass, momentum, and energy. The 
velocities of these particles are drawn from the lo- 
cal Maxwellian velocity distribution function at the 
end of each time step At, thus ensuring the time 
evolution of p , u, and T by the Euler equations. 
Obviously, like the PIC method of Evans and Har- 
low (ref. 6 ), this method requires very long compu- 
tation times. In addition, the Eulerian field vari- 
ables are subject to statistical scatter. As Rich and 
Blackman (ref. 7) modified the PIC method and ob- 
tained the FLIC method for considerably improving 
the efficiency, Pullin (ref. 1) also investigated FLIC- 
like methods. In this variant, cells interact with their 
immediate neighbors through exchange of mass, mo- 
mentum, and energy calculated with the Maxwellian 
distribution. Because the functional form for / is ex- 
plicitly known, it is possible to calculate these fluxes 
by using analytical formulas involving error func- 
tions. It is expected that this FLIC variant will re- 
quire less computational time and the unknowns will 
be free of statistical fluctuations. However, as ob- 
served by Deshpande and Raul (ref. 3), this method 
is only first-order accurate in both space and time 
and has the same stability restriction on At as other 
explicit finite-difference methods. The primary rea- 
son for the restriction on At is that this method per- 
mits fluid movement to immediate neighbors only. 
The true physical situation, however, is that the par- 
ticles can move from a cell to any other cell. This is 
particularly important in high-speed flows. It would 
be highly worthwhile to develop a new FLIC variant 
of the KNM to allow multicell interaction, that is, 
interaction between a cell and its first neighbors, its 
second neighbors, and so forth. Such a method can 
be expected to possess a much less restrictive stabil- 
ity limit on the time step. 

There is one more reason why multicell inter- 
action is desirable. It is well-known that almost 
all explicit numerical methods are very efficient in 
eliminating small-scale errors (i.e., errors with wave- 


lengths of the order of mesh size) but are very slow 
in eliminating large-scale errors. This is the pri- 
mary reason behind the slow convergence to the 
steady state. A FLIC method involving multicell 
interaction is expected to converge rapidly to the 
steady state. The kinetic theory based fluid-in-cell 
(KTFLIC) method of Deshpande and Raul (ref. 3) 
allows multicell interaction. The transport of mass, 
momentum, and energy from a cell to any other cell 
can be obtained in closed form by integrating an ap- 
propriate flux relation involving the local Maxwellian 
distribution. The KNM of Reitz (ref. 2) also involves 
interaction between a mesh point and any other mesh 
point. This method is very similar, in principle, to 
the KTFLIC method and it involves integration over 
velocity only. The KTFLIC, on the other hand, re- 
quires integration over both velocity and space. The 
KNM of Reitz is more efficient than the methods of 
references 1 and 3. 

The KNM’s of Reitz (ref. 2) and of Deshpande 
and Raul (ref. 3) are explicit and unconditionally sta- 
ble. Preliminary computations on a one-dimensional 
shock-tube problem (refs. 2 and 3) show that the nu- 
merical solution does not exhibit pre- and post-shock 
oscillations. There is one major disadvantage, how- 
ever, in the methods of references 1, 2, and 3, namely, 
the magnitude of the numerical diffusion is propor- 
tional to the time step. This restricts the use of a 
large time step in spite of the unconditional stabil- 
ity. Further, with the true physical diffusion being 
proportional to the mean collision time, any exten- 
sion of the method to the viscous case will imply that 
unusually small values of the time step are required. 

The present paper studies the KNM of refer- 
ences 2 and 3 from a theoretical point of view and 
demonstrates that the method is upwinding, satis- 
fies the entropy condition, and is total variation non- 
increasing (TVNI). The entropy condition is satis- 
fied because of the existence of a convex //-function 
defined by Boltzmann in establishing the famous H- 
theorem (ref. 8 ). Further, the question of numerical 
diffusion is studied in considerable detail. It is shown 
that if the basic premise of the distribution func- 
tion is replaced by the Chapman-Enskog distribution 
(with a suitably defined stress tensor r and heat flux 
vector q), then the KNM can be made second-order 
accurate in space and time. The numerical diffu- 
sion in this new variant is of the order of A/ 2 . The 
use of the Chapman-Enskog theory in developing the 
second-order accurate KNM for the solution of invis- 
cid flows is very surprising. The Chapman-Enskog 
theory is always associated with the Navier-Stokes 
equations. In fact, Pullin (ref. 1) has referred to the 
use of the Chapman-Enskog theory for the treatment 
of viscous and heat-conduction effects. Reitz (ref. 2) 
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has used the Chapman-Enskog distribution to correct 
for diffusional effects. The present paper, on the 
other hand, uses the Chapman-Enskog theory to de- 
velop a more accurate KNM for the purpose of solv- 
ing inviscid flow problems. The Chapman-Enskog 
terms are antidiffusive in nature and are a pertur- 
bation over the Maxwellian distribution so that the 
large viscosity of the first-order KNM is cancelled by 
them. Finally, the KNM is extended to the case of 
steady, inviscid supersonic flow. Because of the hy- 
perbolicity of the equations in this case, the KNM 
becomes an Euler space-marching method. 

Symbols and Abbreviations 


a 

speed of sound 

B 

3-7 

4 ( 7 - 1 ) 

C 

normalized peculiar velocity 

c 

peculiar velocity 

c 

_ V2. 

v-i 

D f 

degrees of freedom of a 
molecule 

e 

specific internal energy 

F 

Maxwellian distribution 

F 

contracted local Maxwellian 
distribution 

f 

velocity distribution function 

7,7 

distributions related to / 

7 

contracted velocity distribu- 
tion function 

H 

Boltzmann //-function 

H v 

flux of H 

I 

internal-energy variable due 
to nontranslational degrees of 
freedom 

Io 

internal energy due to non- 
translational degrees of 
freedom 

J 

Boltzmann collision term; 
number of mesh points 

KNM 

kinetic numerical method 

Mi 

Mach number 

n 

time level 


P, 

Chapman-Enskog polynomial 
corresponding to q 

Pr 

Chapman-Enskog polynomial 
corresponding to r 

P 

pressure 

q 

heat flux vector 

R 

gas constant per unit mass 

*D 

ratio of backward to forward 
difference of / 

T 

temperature 

TV 

total variation 

TVD 

total variation diminishing 

TVNI 

total variation nonincreasing 

t 

time coordinate 

t A 

time up to which solution is 
advanced 

U 

conserved variable 

u 

fluid velocity vector 

Ul,U 2 

components of fluid velocity in 
directions 1 and 2 

V 

molecular velocity vector 

Vi,V2,v S 

components of molecular 
velocity in directions 1, 2, 
and 3 

x u x 2 

axes in directions 1 and 2 

X 

position vector 

Xl,X 2 

coordinates in directions 1 
and 2 

X \,L, x l,U 

lower and upper bounds on x\ 

X 2,L, X 2,U 

lower and upper bounds on x 2 

a 

time constant characteristic of 
collisions 

0 

inverse of 2 RT 

1 

ratio of specific heats 

At 

increment in t 

Ax 

increment in x 

6 

angle of deflection 

P 

mass density 

T 

viscous stress tensor 

$ 

multiplier of Maxwellian 
distribution in Chapman- 
Enskog theory 


<t> 

xl) 


general function 

moment function (see eq. (80)) 

collisional invariants 


where I Q = ^ RT for Dj = 1. A solution 

scheme for equation (3) can be constructed by split- 
ting the equation as follows: 


Subscripts: 

CE Chapman-Enskog 

HP high-pressure side of shock 

tube 


df _ _ v df 
dt V dx 


( 6 ) 




(7) 


i,j mesh points 

LP low-pressure side of shock tube 

Superscripts: 

n + 1 intermediate level between n 

and n + 1 


One prime indicates the first derivative and two 
primes indicate the second derivative of a quantity. 


Kinetic Numerical Method 

For studying the basic aspects of the kinetic 
numerical method (KNM), only one-dimensional flow 
is considered now. The Boltzmann equation (1) 
reduces to 

| + (3) 

The velocity distribution / = f(t,x,v,I) is now a 
function of four variables. The moments p , it, e, r, 
and q are defined by 




oo 


p = p(t,x) = 

/ dv 

dl (/) 


J — OO j 

3 

roc 

poo 


pu = / dv 

/ dl (vf) 

J — oo 

/o 


r J 

[°° J 

v 2 \ 

pe = j dv 

\ dI 

I + T / 

J — oo 

Jo V 

2 / 

rOO 

poo 


p — r = I dv 

/ dl(v 

-u) 2 / 

J —oo 

Jo 


poo 

r oo 



q = - 1 dv j dl 

J ! (v-u ) 2 

(v - u)f 

J —oo 

Jo 

2 



(4) 


In the Euler limit, 


With the distribution function at time level n de- 
noted by f 1l (x,v,I), the solution of equation (6) is 
given by 


f n+1 (x, v,I) — f n (x — v At, v, I) (8) 

Note that equation (8) is an exact solution of equa- 
tion (6), which is a collisionless Boltzmann equa- 
tion. In the Euler limit, J(f,f) = 0. The vanishing 
of J is due to a very large number of intermolecu- 
lar collisions. Because J = 0 if and only if / is a 
Maxwellian distribution, it follows that the solution 
of equation (7) in the Euler limit is 


F n+l {x,v) = F(/ n+1 ) (9) 

where F ^/ n+1 j is the local Maxwellian distribution 
with p, it, and e the same as those of / n+1 , that is, 


// 

The method thus consists of two phases, the convec- 
tion phase and the collision phase. In the convec- 
tion phase, the fluid particles located at x — v At are 
moved to x, whereas in the collision phase the re- 
sultant distribution function after convection imme- 
diately relaxes to the local Maxwellian distribution. 
In view of the conservation of mass, momentum, and 
energy during the collision phase, the first three mo- 
ments of F n+l are the same as those of f n+i . There- 
fore, the KNM can be written as 


1 

V 2 
r i JC 

. ' IT J 


(F B+1 - f n+1 )dvdl = 0 (10) 


J(fJ) = 0 

and the distribution / is then the Maxwellian 
distribution: 


f=-F = 


P 

I 0 {2ttRT) 1 / 2 


exp 


(v — u ) 2 
2 RT 



p n+1 = 
3 

poo 

j dv 

poo 

/ dl F (tn , Xj - 

vAtj) 


J —oo 

*0 



poo 

poo 


(p«)” +1 = 

j dv 

/ dl vF (tn , Xj 

— v At, 


J —oo 

1 0 


poo 

r°° / ..2 \ 


(pe)" +1 = 

I dv 

/ — oo 

/„ 

*3 

f 

H 

1 


( 11 ) 
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where the standard notation pd = p (t n , Xj) has been 
used. The integration with respect to the variable / 
yields 


rOO 

p n+1 = dvF (t n , Xj -vAt) 


J — oo 



rOO 



(pu ) n+1 = / dvvF (tn,Xj ~ v At) 


J — oo 



/•OO 

r ..2 1 


(„e )" +1 = / dv 

+ 

< 

a 

1 

''-I 

H 

j: 


J — OO 

L J 


X F ( tn,Xj - v At) 



(12) 


where F is the contracted local Maxwellian distribu- 
tion defined by 


F = p 



~P{v — w) 2 


(13) 


where (3 = Equations (12) are basic relations of 
the KNM. The field variables p, u, and T are stored 
at mesh points only, and therefore F (xj — v At) 
has to be determined by some kind of interpolation 
(unless xj — v At is chosen to coincide with another 
mesh point). One choice for interpolation is 

F (x 3 - v At) = F t (l -0 + F i+l £ (14) 

where Xj — v At = + £ At for 0 < £ < 1. The 

relationship between x t , Xj, and £ is shown in the 
sketch below. 


W- 


-v At- 


r I • 


i' 

-*J 

I 


i+1 


N — 

V,. 


M 

I 

— 4 - 

x. , x. 
J-l 3 


% Ax 


x j+l 


It is possible that the point Xj — v At is outside 
the computational domain, that is, the fluid particle 
has crossed the domain boundaries during convec- 
tion. Obviously, a boundary strategy is then required 
to find F(xj — v At) so that physically meaningful 
desired boundary conditions are satisfied. A discus- 
sion of how to specify boundary conditions is given 
later. 

The interpolation formula (eq. (14)) is only first- 
order accurate in space and equation (8) is only first- 
order accurate in time. Hence, the KNM (eqs. (12)) 
for the solution of the Euler equations is first-order 
accurate in space and time. It is shown subse- 
quently that to achieve second-order accuracy it is 


not enough to use a second-order accurate interpola- 
tion formula. 

A few important properties of equations (12) are 
worth noting. The implementation of equations (12) 
involves only quadrature, and no numerical differen- 
tiation is required. Because of the long-range inter- 
action among mesh points the method, despite be- 
ing explicit, is expected to be unconditionally sta- 
ble. Indeed, Reitz (ref. 2) and Harten, Lax, and van 
Leer (ref. 9) have shown that methods of the type 
of equations (12) are unconditionally stable. Also, 
the method possesses a high degree of vectorization 
because it is explicit and dominated by arithmetic 
operations rather than logic. The vector code is, in 
fact, about 8 times faster than the scalar code for a 
one-dimensional shock-tube problem with 500 mesh 
points. 


Properties of the KNM 

At this point, it would be interesting to investi- 
gate whether the KNM has other important prop- 
erties like upwinding, entropy condition satisfaction, 
and TVNI. Of late, development of schemes possess- 
ing these properties is considered highly desirable be- 
cause of many advantages. For example, an upwind 
method is robust and has a lot of physical appeal 
(ref. 9). A method satisfying the entropy condition 
prevents formation of expansion shocks. Further, if 
a method is TVNI then the solution is free from pre- 
and post-shock oscillations. In the case of the KNM, 
the upwinding nature is a direct consequence of equa- 
tion (8), which states that the solution at any mesh 
point is influenced by the data upwind of that point. 
As mentioned in the Introduction, an //-function can 
be defined which satisfies the entropy condition. The 
demonstration of this capability proceeds as follows. 

The Boltzmann //-theorem has been described in 
the kinetic theory of gases as the bridge connecting 
the equilibrium thermodynamics with the nonequi- 
librium statistical mechanics. Briefly stated, it says 
that the //-function defined by 


H = III f]nfdv (15) 


monotonically decreases with time as a homogeneous 
gas in statistical nonequilibrium evolves to equilib- 
rium. In the case of spatial inhomogeneity, the theo- 
rem states that 


dH 

dt 


^ dH { n 

+ Eis 7 s0 
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where H{ is the H- flux defined by 

H { = fjj Viflnf dv i dv-i dv 3 

For the one-dimensional case we therefore define 
5-37 


-in 

“II" 


H 

Hy = 


F\nF + 


F\nF + 



(16) 


2(7-1) 


Note that an additional term involving F\n(3 is 
required to account for nontranslational degrees of 
freedom. (See ref. 8.) From equations (16) we obtain 


dH 

dt 


dF dF\ . 

— h v — — dv dl 

dt dx ) 


d -£dh'° F >( 

■JHb-h) 


5-37 


2(7 - 1 ) 

x (Fin ())dv dl 


(17) 


Now, 


/« (f + ”§r ) ' ■ • /“ [I M + 4 ('“'■)] d " 

(18) 

where F is the contracted local Maxwellian distribu- 
tion defined in equation (13). Substituting for I 0 in 
terms of (3, we get 


//*(£♦•£)“-/ 


3-7 
4(7 - 1 ) 


d F 
+ v — ( — 
dx \ 0 


3-7 f 

*(Tf - 1) J 


4(7 
- F 
3-7 


d_ ( F 
di [p 

dv 


dF dF 

(- v 

dt dx 


( i lD 0 +v ic ]nl3 )] dv 


I 


(1 + In P) 


4(7 - 1) 

'dF dF , , 

X I h v I dv 

dt dx 


3-7 
4(7 - 1) 


/ [- 
J iat 


(Fin p) 


+ v — (Fin P) 

OX 


(19) 


Noting that the equation of continuity is 


we find equation (19) yields 


II 




v ) dv dl = P — -H 

dx j 4(7 — 1 ) 


/ 




+ v — (Fln/3) dv 
dx 


d_ 
dx 

3-7 
4(7-1) 

J t (FlnP) 


IS 


+ V — (Fln/3) 
dx 


dvdl (20) 


With equation (20), equation (17) gives 


dH dH, 
dt dx 


?- in 


1 + In F — 


10 — 67 
3-7 


an 


\ dt 

If we observe that 


l dF , 

X ( + v — — ) dt) d/ 


dF\ 

’dl) ' 


( 21 ) 


IO-67 / 3 4(7 - 1 ) 

InF- lpi={ lnp + -ln/3 + ln— £ L- - pu‘ 

3 — 7 y 2 7 r 1 /^(3 — 7) 


+ 20uv 


- m ( , + t) 


( 22 ) 


is the linear combination of collisional invariants 1, v, 
2 

and I + \ and further that the Euler equations of 
motion can be cast as 


II 


/dF dF\ 

,p {-^ +v w dvdI=0 


ip = 1, v, I + 


equation (21) reduces to the entropy conservation 
dH . dH v 


(23) 


dt dx 


= 0 


(24) 


Thus, every smooth solution of the Euler equations 
satisfies the entropy conservation equation (24). We 
will now demonstrate that H decreases during the 
collision phase in conformity with the H- theorem. 

During the collision phase the dis tribu tion func- 
tion suddenly relaxes to F n+1 from / n+1 = f n (x — 
uAf). Hence, the jump in H during the collision 
phase is given by 


H n+1 


H n+1 = j j ^/ n+1 ln/ n+1 -F n+1 lnF n+1 ) dvdl 

+ k jj (/"+ T ln/3"+ T 
- F n+1 ln/3 n+1 )di;d/ (25) 


6 



where k = ^-l) and ^ n+1 = ^ n ( x _ vAt). 
Considering the first term on the right-hand side of 
equation (25), we have 


f [ (/ n+1 ln/ n+1 - F n+1 lnF n+1 'j dvdl 

- IS ^ 


pn+l 


dv dl 


+ 


3-7 


if ( fn+1 ~ Fn+1 ) ]nFH+1 dvdI 

= if f^lnt^Ldvd!* 10 -^ 

X II 0 n+1 I ( 

- II r* 


rfi+1 ipn+1 


dvdl 


fTl+1 

In - — r-rdvdl 


pn+l 


IO-67 

3-7 


j 0 n+1 (/o +1 / n+1 - /? +1 F n+1 ) dv 


_ _ ( 26 ) 

where f n+1 and F n+1 are contracted distributions. 
Through substitution for /„, the right-hand side of 
equation (26) becomes 


because of equation (10). Equations (25), (27), and 
(28) finally yield 



The second term on the right-hand side of equa- 
tion (29) is positive because of the inequality 

x— 1 > lnx . (x > 0) (30) 

The first term can be shown to be positive by pro- 
ceeding as follows. We have 

// * {Jr* (7^ - : ‘) 

<jj (f" +1 - / n+1 ) dvdl 
< 0 


/ n+1 -F n+1 j dv 


IK"*'' n 7 ^ d ' ,J,+ K& 

Noting that 

p n+l = JJ pn+l dvdJ _ J p>n+l dv = J J f n + 1 dvdl 


=! 


pn+l 


dv 


we obtain 


//( 


/ n+1 ln/” +1 - F n+1 In F n+1 \ dvdl 
dv dl 


If 


f n+l In fr>+l 


pn+l 


^II r+ '&- I ) d ^ (27) 

For the second term on the right-hand side of equa- 
tion (25), we have 


■IK 


f n+1 ln0 n+1 - F n+l In 0 n+1 J dvdl 
+ k j j ln/3 n+1 ^/ n+1 - F n+1 j dvdl 

-'ll ^ 


■n+T,_ P 1l+1 
0 n + 1 


dvdl 


(28) 


because of equation (10). Hence, 

(31) 

According to Kullback (ref. 10), the right-hand side 
of equation (31) is the information theoretic dis- 
tance between two distributions f n+ 1 and F n+ 1 . 
This distance is always positive and vanishes only if 
fn+l _ pn+ l Combining equations (29) and (31), 
we determine that the //-f unction decreases in the 
collision phase from H n+1 to H n+l . This sudden 
reduction in the //-function is due to the instanta- 
neous relaxation of the velocity distribution function 
from F n {x — v At) to F n+l , and is further due to 
the convexity property of the //-function. 

The validity of the TVNI property for the KNM 
can be established as follows. Let the distribution 
function at a mesh point j at the nth time step be 
denoted by F' 1 . Then, 

/p = (i - 0*7* + ^+1 (0 < z < 1) 

which is a special case of the slightly more general 
interpolation scheme 

P = (1 - q) Ff + c/f +1 (0 < < 1) (32) 
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Here i is the mesh point such that Xj — v At lies 
between x l and Xj +1 . The variable ej can depend on 
the particle velocity v. The total variation of / is 
defined by 

TV(D = £ |/f +1 - ff | = £ |*7 +1 - (33) 

all j allj 


During the collision phase, the density, fluid velocity, 
and temperature do not change (eq. (10)). Hence, 


TV (p 


n+l\ _ 


TV 

TV 


(»“ +1 ) 


= TV 


TV 


(/■ + ‘) ) 
(” n+1 ) 


_ 'py ^yn+1^ 


(34) 


Therefore, it is sufficient to deal with the change of 
total variation during the convection phase only. The 
total variation of / n+1 is given by 

tv (r+i) = £ - /” +1 1 < £ i (i 

allj all i 

+ £‘ l+J i^ 2 -^ii 

all i 

Evidently, 

E'i+ii^+z - Ail = E f iic n +i - F ? 

all i all i 


)\n+i~ F i\ 


(35) 


and hence the inequality (35) becomes 


TV ( f n+i ) < £ I F? +l - F ? I < TV (/ 


all i 


linear, hyperbolic equation with the vector conser- 
vation law by use of the fact that the law is a mo- 
ment of the equation. From the above analysis the 
TVNI property of the KNM hinges upon the interpo- 
lation strategy adopted in calculating F(x - v At). 
This fact will be used later while developing a second- 
order accurate TVD KNM. 


A Second-Order Accurate KNM 


It has been mentioned in the Introduction that 
the KNM’s of Pullin (ref. 1), Reitz (ref. 2), and 
Deshpande and Raul (ref. 3) suffer from a major 
disadvantage in that they have numerical diffusion 
proportional to the time step. From the physical 
point of view, such a result is only to be expected 
because the fluid particles in the KNM are allowed 
to move over time step At before they undergo 
collisions. The distance traveled between collisions 
is thus proportional to At. From the kinetic theory 
it then follows that the mean free path, and hence the 
viscosity, will go like Af. This is a very large amount 
of viscosity (as the results shown later will verify). 
Therefore, a modification in the KNM developed 
earlier is required that will ensure that the method 
has higher order numerical viscosity. The analysis 
of this section shows that it is possible to achieve 
this aim by the use of the Chapman-Enskog theory. 
It becomes clear later that the development of the 
second-order accurate KNM is intimately connected 
with obtaining higher order numerical viscosity. 

For one-dimensional compressible inviscid flow, 
the governing Euler equations are 


with the use of equation (33). Thus, for every value 
of v and I, the total variation of the velocity distribu- 
tion function is nonincreasing during the convection 
phase. There are some advantages to considering the 
TV of / instead of the TV of p, u, and T. The vec- 
tor conservation laws of the Euler type are moments 
of a single scalar equation for /. This equation be- 
ing linear and hyperbolic admits the exact solution 
f(x — v At), which is TV-preserving. It is therefore 
quite consistent to demand that the numerical so- 
lution for / be TV-preserving. Further, the TVNI 
condition on / leads to a wiggle-free solution for the 
vector conservation equation. On the other hand, the 
TVNI condition on each of the vector components 
p, pu , and pe is, in general, of dubious physical va- 
lidity even though it gives wiggle-free solutions. The 
present approach connects the TVNI property of the 


+ = 0 

|w + ^(P + f» 2 ) = ° • 

|w + ^(^ + p») = o 


(36) 


As Df = 1, equations (4) and (5) yield 


r RT u 2 

e = Io + ^r + t - 


RT u 2 
7- 1 + Y 


(37) 


The second-order accurate Taylor expansions for 
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p n+1 , (pu) n+1 , and ( pe) n+l 


are 


F dvd.1 = p 


' n+1 =' B+ (ii) At+ (^J + ° (At3) 

pu )" +1 = (pu)" + (^pu) At 

+ T i+0(Al3) 


(pe )" +1 =(pe)”+ (£pe)’ 


+ (S'*)"1 !+0(a ‘ 3) 


vFdvd.1 = pu 


v 2 F dvdl = p + pu 2 


F dvdl — pe 


v 3 F dv dl = 3 pu + pu 3 


v 2 FdvdI = -p- -pu 2 


2(T-1)' 


T / P 2 \ , P« 4 


+ ^ri7; + » 


These expansions contain the first- and second-order 
time derivatives of p , pu, and pe. The first-order 
time derivatives can be replaced in terms of the first- 
order space derivatives by equations (36). Replac- 
ing the second-order time derivatives in terms of 
the space derivatives requires detailed manipulations, 
which are shown in appendix B. With the results of 
appendix B, the Taylor expansions (38) become 


p"+l =p"-At£(pu)" + ^|i(p + pu 2 )" + 0(At 3 ) (39) 


Through use of the moment equations (42), the 
Taylor expansions (39), (40), and (41) can be cast 


= \ pu \ — At- 


a F f 
2 J 


dv dl vF n v 


dvd!v 2 F n v 2 


^—(3 -l)P& 

¥(3-a)^ + ^ p£A E p). 


+ 0(At 3 ) 


(pu ) n+1 = (pu)" - a t^-( P +pu 2 r + + p« 3 ) n 

+ °^ ( 4 °) 


From use of the second-order accurate Taylor 
expansion 


dF\ n v 2 A t 2 ( d 2 F 


(3 F \ 71 

— J + 


(pe)"+ X = (pe)" - At — (peu + pu)" 


At 2 d 2 5~f — 3 9 7 p 2 pu^ 


At 2 a 


9 h, , 9u 7 9 ( p\l n 


+ 0(At 3 ) 


and with the definitions 


At, du 

r = - T (3-9)p^ 

Af 7 d ( p 
q ~ 2~ 7 - l P dir \p 


the Taylor expansions (43) take on the simple form 


= JJ dvdIF n (x — v At) v 2 

U + V. 


tu + q 


r +0{At 3 ) (45) 


Notice that the local Maxwellian distribution F sat- 
isfies the moment equations 
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Equation (45) reveals a very important feature of 
the KNM. The right-hand side of equation (45) con- 
sists of two terms; the first is an upwind term be- 
cause of the presence of F n (x - vAt), and the sec- 
ond is an antidiffusive term containing r and q. The 
antidiffusive nature results because r and q defined 
by equation (45) are analogous to negative viscous 
stress and negative heat flux. To achieve second- 
order accuracy in both space and time, it is essen- 
tial to consider both the upwind term and the an- 
tidiffusive term. Further, F(x — v At) appearing in 
the upwind term must be evaluated for second-order 
accuracy in space by using a suitable interpolation 
scheme. Note that the antidiffusive term is 0(At 2 ), 
and therefore second-order accuracy in space is at- 
tained with the use of the upwind term alone if At is 
chosen such that O(Af) = 0(^Az 3 / 2 ^. For Courant 
numbers of the order of unity, the antidiffusive term 
is 0( Ax 2 ). Thus, any version of KNM in which only 
the upwind term is considered and the Courant num- 
ber 0(1) is used cannot achieve second-order accu- 
racy regardless of the accuracy of the interpolation 
scheme. This fact has not been noticed in the earlier 
work. For example, Reitz (ref. 2) ignores the anti- 
diffusive term and uses a second-order accurate in- 
terpolation formula for calculating F(x — v At), 
presumably to achieve second-order accuracy. As the 
Courant number chosen by him is 0(1), it is clear 
that his calculations are subject to the truncation 
error 0( Ax 2 ) and thus are first-order accurate. 

The preceding analysis shows that a modification 
in the ansatz is necessary for constructing a second- 
order accurate method. This modified ansatz f 0 — 
fo(t,x,v, I) must satisfy the constraints 


If 


1 

v 

I + % 


-// 


2 J 

1 

V 

7+^ 


L d_ 

dx 


fo(t, x - v At, v, I)dv dl 

F(t, x — v At, v, I)dv dl 

+ 0(Af 3 ) (46) 


0 

t At 
[{ru + q)At\ 



1 



v 2 (f o -F)dvdI = 0(At) 


(49) 


The conditions in equation (47) merely state that 
the density, velocity, and energy for f Q and F must 
be identical within the truncation error. They are 
therefore conservation conditions. With the defini- 
tions in equations (4) for viscous stress and heat flux, 
the conditions in equation (48) can be equivalently 
cast as 

JJ ( v — u) 2 f 0 dvdl — p — t (50) 



/ + 


( v — u ) 2 
2 


(u - u)f 0 = -q 


(51) 


Equations (50) and (51) require that the distribu- 
tion function f a must have nonvanishing viscous 
stress and heat flux. They therefore suggest the use 
of the Chapman-Enskog distribution instead of the 
Maxwellian distribution. This is further supported 
by the observation that equation (43), which is the 
basis of equation (46), is obtained by replacing the 
time derivatives of field variables with their space 
derivatives using Euler equations. Such a replace- 
ment is very characteristic of the Chapman-Enskog 
analysis. 

Consider the model Boltzmann equation 




(52) 


where a is a time constant characteristic of collisions. 
The parameter a may depend on space but is inde- 
pendent of v and /. For large values of a it is advan- 
tageous to write equation (52) in the form 


/ = 





(53) 


The Chapman-Enskog distribution is given by 


These constraints are satisfied if 

(fo- F)dvdl = o(At 3 ) (47) 


1 

v 

U + ? 


II 

//Lb 


F) dv dl = 


0 
T 

(ru + q) J 


+ 0 (At 2 ) 

(48) 


, „ 1 (dF dF\ 

fcE = F ~a (W +V te) 

With the results of appendix C, the above distribu- 
tion can be written as 


/CE = F 




(54) 
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where 


r CE 

QCE 


a ox 
7 d 
0(7 — l)^ dx 



> 


(55) 


Pr = 

P„ = 


3t — 5 
2(3 - 7 ) 

2(7-1) 

7 


4(7 -l) 2 * C 2 

3-7 2PT 


5C ^3 

+ C 3 

2 


4 (7 ~ 1) c I 
3-7 2PT 


(56) 

(57) 


(2RT) 1 / 2 

Equations (55) and (44) immediately fix a at a = 
— ~Kt' The model equation (52) then becomes 


dt dx A t y J) 


The characteristic collision time for this equation is 
half the time step At, and this result merely confirms 
the earlier argument about the numerical diffusion 
being proportional to At. The sign of a is negative 
because the numerical diffusion inherent in the first- 
order KNM has to be cancelled to develop a second- 
order KNM. Thus, the Chapman-Enskog terms in 
equation (54) are antidiffusive. That fcE defined 
by equation (54) is the required distribution can be 
clearly seen from the following relations: 


Several important features of the above second- 
order accurate KNM are worth noting. Equa- 
tion (58), containing the Chapman-Enskog distribu- 
tion, has been derived from equation (45). As stated 
previously, the right-hand side of equation (45) con- 
tains two terms. The first term, which is a mo- 
ment of F n (x — vAt), is upwind in character. The 
second term cannot be expressed as a moment of 
F n (x — v At) and is antidiffusive. The antidiffusive 
term may be absorbed in the upwind term only if 
the distribution function is not Maxwellian. Equa- 
tion (58) is an upwind version of a second-order ac- 
curate solution in which the perturbed Maxwellian 
distribution or, equivalently, the Chapman-Enskog 
distribution is employed. Obviously, tqe and qcE of 
the Chapman-Enskog distribution /q E are antidiffu- 
sive. Furthermore, because of the presence of tqe 
and qcEi the Chapman-Enskog distribution not only 
depends on local values of field variables but also 
depends on their neighboring values as well. The 
support of the Chapman-Enskog distribution is thus 
more than that of the local Maxwellian distribution 
used in first-order accurate KNM. 

The Chapman-Enskog distribution function /qe 
can assume physically meaningless negative values 
for sufficiently large values of velocity v. From 
equation (54) it is obvious that /q E can become 
negative if 


<7CE 

p(2RT) 1 / 2 


= 0 ( 1 ) 


(59) 



(/ce — F)dvdl = 0 



v ~ u) 2 f C E dvdl = p - r CE 



(v-u)f C E dvdl = - q CE 


(fCE~F)dvdI = 0(At) 


These moment equations can be easily derived by 
integrating /q E . They are identical with the condi- 
tions (47), (48), and (49). It therefore follows that a 
second-order accurate solution in time is given by 


p 

pu 


-i n+1 


Lpe J 



+ O (At 3 ) 


v At, v, I)dv dl 


(58) 


that is, if v > O (^At 1 / 3 ^. As Narasimha (ref. 11) 
has noted, the Chapman-Enskog distribution is only 
an inner solution in u-space and is not valid for 
relatively high velocities. It is possible to modify 
the Chapman-Enskog distribution along the lines of 
Narasimha’s analysis (ref. 11). 

Another important point about equation (58) is 
that fcE( x j ~ v At) needs to be evaluated at vari- 
ous values of v. Hence, as noted before, some kind of 
interpolation scheme is required. This scheme must 
be second-order accurate in space, should yield non- 
negative interpolated values, and should satisfy the 
TVNI sufficiency condition of equation (32). These 
questions are addressed in a subsequent section. 

Interpolation Scheme 

The central problem in devising an interpolation 
scheme is to evaluate f(xj — v At), given the values 
of / at mesh points, so that the conditions set 
forth previously are satisfied. Let the mesh point i 
corresponding to point j be such that 
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Xj — v At = xi + £ Ax (0 < £ < 1) (60) 


The function f(xj — v At) = f(x j + £ Ax) then 
depends on the neighboring mesh points i ± 1. In 
fact, the Taylor expansion yields 

c t2 

f ( x i + i Ax) = fi + - (/i+l ~ fi-l) + — (fi+1 ~ 2 fi + fi- 1) 


Unfortunately, the above expansion does not auto- 
matically ensure positivity of f(x ? + £ Ax) even if fi 
and fi ± i are assumed to be positive. This is particu- 
larly true for calculations near shocks. To devise the 
desired interpolation formula, first write the above 
Taylor expansion equivalently as 


/Ui + $ Ax) 


Now, 


( 1_ 0/» + £/t+l - — 9 — (/i+l “ 2 /i + fi-l) 

(61) 


^ (A+l-2/< + /i-i) 


e(i-o 2 


(/i+l - 2 /i + / i— 1 ) 


e 2 (i-o 


(/i+l “ 2 /i + /i — 1 ) 


= 2e(1 - 0 2A + l- 2 /i+^ l/ 
v ’ k+i- 2 k + k-i 1 

+ 2e 2 (i - ^ + o(Aa;3) 

= 2£(1 - £)fc [(1 - 0/t + e/i+i] + 0(Ax 3 ) 


where 

, _ /i+l ~ 2 /i + fi-l 
1 /i+l + 2/i + /i-1 

The function fa has the property \fa\ < 1 if f l and 
/ ( +] are positive. In terms of fa the interpolation 
formula (eq. (61)) becomes 


with the constraints 0<ej<l,0<q<l, and 
Cj + e'i < 1. With the method of Chakravarthy and 
Osher (ref. 12) limiting the contribution of the second 
difference, it is possible to devise an interpolation 
scheme that satisfies the TVNI sufficiency condition. 

The second-order accurate Taylor expansion can 
be written as 


/ (x t + £ Ax) = fi + | (/ i+1 -fi + fi- /i_i) 


where 


+ 

= f% + 
+ 


Y (/i + l ~ fi + fi-l ~ fi) 
5(1 + 0 (/*., fi) 


2 

S(1 - 0 


{ft -fi-l) 


«1 + 0 


— fi + (/i+l ~ /i) 

+ o 


Backward difference /j — / t _i 

D Forward difference fi+i ~ fi 

In smooth regions, 


(63) 


(64) 


r D = 1 - Ax J -j, + 0(Ax 2 ) 

* i 

and thus rjp remains close to unity. In flow regions 
near shocks or contact surfaces, rjj can wildly vary 
and some limiting criterion is required to preserve 
the TVNI condition. The key to satisfying the TVNI 
condition lies in requiring 


0 < — r D < 1 (65) 


Two cases arise, namely, > 0 and tq < 0. If 
we consider the case of r/j > 0, the condition in 
equation (65) is satisfied if 


/(*i + € Ax) = [(i - Ok + €/i+i] [1 - 2£(1 - 04>i] (62) 

Because |2£( 1 — £)£q| < 1 if f, and / t±i are posi- 
tive, it follows that equation (62) yields only positive 
interpolated values as long as the distribution func- 
tion is positive. Unfortunately, equation (62) does 
not satisfy the TVNI sufficiency condition. This is 
because fa can become negative, and in that case, 

i-2£(i-o& = i + 2e(i-ow>i 

It is not possible then to write equation (62) in the 
form 

f{xi + £ Ax) = eJi + e'ifi+i 


vd < 1 + | ( 66 ) 

As 0 < £ < 1, the right-hand side of the above 
inequality has the minimum value of 3. One way 
of satisfying equation (66) is to limit the value of rjp 
to 3 when > 0. For the case when rp < 0, the 
condition in equation (65) is satisfied if 

or, equivalently, if 

\ t d\ < (67) 
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Thus, by limiting the relative values of the forward 
and backward differences and taking rp> = 1 out- 
side these limits, we find the interpolation formula 
in equation (63) yields not only positive values of 
f(x{ + £ Ax) but also satisfies the TVNI condition. 
As mentioned previously, the basic input to equa- 
tion (63) is the set of positive values of /j at all mesh 
points. 

From the analysis in previous sections it has be- 
come clear that the second-order accurate KNM re- 
quires the use of the Chapman-Enskog distribution 
and a second-order accurate interpolation scheme. 
The development of the KNM so far has been re- 
stricted to one-dimensional unsteady, inviscid com- 
pressible flows. The extension of the KNM to two- 
dimensional steady supersonic flow is considered in 
the next section. 


are related to / through the moment equations 


P = 

pu - 


III 

III 


f dv\ dv 2 dl 
vf dv\ dv 2 dl 



f dv\ dv 2 dl 


(71) 


The KNM for equations (68) can be constructed by 
splitting the Boltzmann equation (69) as follows: 


df df 

1- v 2-^— = 0 

OX 1 OX 2 


(72) 


0 (73) 


Extension of the KNM to Two-Dimensional 
Steady Supersonic Flow 

So far, the development of the kinetic numerical 
method (KNM) for the unsteady Euler equations has 
been considered. The steady Euler equations for 
supersonic flow are hyperbolic, and hence a space- 
marching version of the KNM can be developed along 
similar lines. In this section, various steps involved in 
developing the Euler space-marching KNM are given. 

The two-dimensional steady Euler equations are 


d d 

^-(m) + ^(p U2 ) = ° 

a|-(p + ^) + ^(p Ui “ 2) = 0 

d d 

^-(^^2) + — (p + p U |)=° 

(peiti +pui) + ( peu 2 + pu 2 ) = 0 


(68) 


where e = -r-2 — . + + u 2 

wncie e p ( 7 _i)T 2 

Equations (68) can be obtained as the moments 
of the steady Boltzmann equation 


ui 


df_ 

dx\ 


+ V2 


dj_ 

dx 2 


= J(U) 


(69) 


if the velocity distribution function / is Maxwellian: 


F = 


2irRTI 0 


exp 


(vi - Ml) 2 + (V2 - U2F 

2 RT 



where I 0 = ^ J RT. The field variables p, u, and T 


The exact solution of equation (72) is given by 


f (x i + Axi,z 2 ,i>i,H2,7) = f (x 1,12 - c Axx, Aii,v 2 ,/) 

(74) 

where c = In the convection phase, therefore, the 
fluid particles at x\ + Aaq and x 2 come from x\ and 
x 2 — cAx\. The collision phase, governed by equa- 
tion (73), merely says that the distribution function 
/ is a local Maxwellian distribution everywhere. 

Consider a two-dimensional mesh of points (n,j) 
within a rectangular computational domain. The 
index n is along the Xi-axis and the index j is along 
the X 2 -axis. Because x\ is a time-like coordinate in 
the integration of the hyperbolic equation (72), the 
finite-difference solution for the line n + 1 can be 
obtained from the data for the line n. We introduce 
the following notation: 

f n = f{n Axi, x 2 , vi, v 2 , I ) = f{n Axi, x 2 ) 


The solution (74) can then be written as 


fj +1 = f(nAxi,x 2j -cAx 1 ) (75) 


(p«i) 


= Iff vif (n Axi, x 2 j - c Axj) dvi dv 2 dl (76) 
(p + P“i)" + = III V 1 f ( n A xi , x 2 j — c Axj ) dvi dv 2 dl 


(77) 


{pu iu 2 ) j n+1 = I II Viv 2 f (n ^xi,x 2 j - c An) dvi dv 2 dl 

(78) 


[peu\ + pui 


>r=UJ 


n |/+y 


x / (n Ax i , x 2 j — c Axi ) dvi dv 2 dl 

(79) 
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Note that because of equation (73) the velocity dis- 
tribution / is the Maxwellian distribution in equa- 
tion (70). The moment function ' 4 ' is defined as 


Thus, P CE = g^-(lnF) + c^(lnF). Replacing RT 
in equation (70) by we get 


<4 = 



(80) 


and the inner product is defined by 


F = 


p 3 (7 ~ i) 

27 rp^ (2 — 7 ) 


X exp 


p(v\ - u\) 2 + p(y 2 ~ u>2) 
2 P 


pl 1-1 

P 2-7 


ipF(v,I)dv 1 dv 2 dl 


with suitable limits for iq, iq, and /. Equations (76) 
to (79) can then be written in the following com- 
pressed notation: 


which gives 


In F = 3 In p — 2 In p — 


pi ~t - l __ P (vi - «l) 2 + P (^2 ~ ^2) 2 
P 2-7 2p 


+ In 


7-1 
2ir(2 - i) 


= (9,f(x i n , -cAzq, /) ^ (81) 

which is obtained by taking the ^-moment of equa- 
tion (75). The KNM space-marching scheme for 
the Euler equations (68) is therefore simply the \4- 
moment of the exact solution of the collisionless 
Boltzmann equation (72). The KNM space-marching 
scheme in equation (81) is first-order accurate in x\ 
and X 2 - 

For the purpose of developing a second-order 
accurate Euler space-marching scheme, we first note 
that 

/ dF \ 

(9,F (x\ + A11, 12)) = (W,F) + Axi 



The ^-moment of the collisionless Boltzmann equa- 
tion (72) yields 



Combining equations (82) and (83) gives 


(®, F (x! + Ax!,X 2)) = (’*>, F) - Ai! 




From equation (84) it is obvious that to develop 
a second-order accurate KNM, the second partial 
derivative of F with respect to x\ must be obtained 
in terms of its derivatives with respect to x 2 . For 
this purpose, we proceed as follows. Let 


dF 

dxi 


.OF 

°d7 2 


= p ceF 


(85) 


Hence, 

\dx! 3x 2 / V <3xi 3*2 / 

( v l ~ “l) 2 + ( v 2 - u 2) 2 


3 
P 

I 7 — 1 
p2 — 7 


2 P 


+ 


p(«i - m) 
p 

( 3^2 j3u2 \ p{v 2 - u 2 ) 
\9xi 3x 2 / p 


/au 1 + .3u 1 \ 
V 3xi 3x 2 / 


V3xi 3x2/ 


_2 p(vi - m) + p(«2 ~ «2) 
P 2 p 2 


pi 7-1 


P 2 2-7 


After a little manipulation we obtain 


p ce = 


= (^ + ~ c ^\ 


\dxi <3x2 7 
1 / 3p 

P \3xi 


In F 


i£. + g i£-') 

3x2 / 


x 3 


(«i - ui) 2 + (n 2 -U2) 2 7-1 I 


2RT 


2-~i RT 


j 3«1 „ 3«i ^ l>i — Uj 1 ( 3«2 , - 3u2 \ v 2 - «2 
+ \ 3x7 ^ C 3*2 / + \3^7 +C 3^J AT 


1 f _3p^ 

P \3x! 


, - 3p \ 

+ C ^J 


I 7 “ 1 + ( V 1 ~ u l) 2 + ( v 2 - u 2) 


RT 2-7 


2RT 


( 86 ) 


The function Pce is a polynomial in c, iq, and /. 
This is evident after replacing t>2 by ciq. The space 
derivatives of p, u\, u 2 , and p with respect to x\ can 
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be replaced by X 2 space derivatives from the Euler 
equations. (See appendix D.) Now 


d 2 F 

9x 2 




dx\ V dx 2 / dx\ 
d 

cPcef) 

9 


9x2 

— r 

9x2 V 


9 f _ 2 9F „ „\ . 9 

c — h 


9x2 




5 2 S-^3?7<Pc Ef) + 


9x^ 9x2 


9xi 


(PceE) (87) 


Notice that 




£<*•**') 

a / dF ~dF 

dxi \ ’ dxi C di2 


= 0 


This is because 


/ dF „dF\ n 
represents the Euler equations. Therefore, 


(88) 

Substitution for in equation (84) finally 


gives 


{^.Ftn + Ax 1 ,x 2 )> =(*,F) - 


Ax? 


Axr Id \ 

■■^r ,2 9^ (PcEF) ) 

+ o(Axf) (89) 

The first three terms on the right-hand side of 
equation (89) can be further simplified with the 
Taylor expansion as follows: 


dF c 2 Ax 2 d 2 F 

9x2 2 


F (x 2 - c Aii) = F - c Axi i- ^-2 + O (Axf ) 


Equation (89) can be alternatively cast in the form 


(*,F(x 1 + Axj , x 2 ) ) 


(*, F(xi,x 2 - c Axj)) 


Axf / 3 \ 

+ O (Axf) (90) 


The inner product of with F (x\ + Axi, X 2 ) 
gives the field variables pu\,p + pu\, puyU 2 , and 
peu\ + pu\ at xi + Axi and X 2 , and hence equa- 
tion (90) constitutes the second-order accurate space- 
marching KNM for the Euler equations. It is in- 
teresting to note that the right-hand side of equa- 
tion (90) contains an upwind term and an antidif- 
fnsive term. The antidiffusive term is absent in the 
first-order accurate KNM (eq. (81)). However, for 
second-order accuracy it is not enough to consider 
E(xi,X 2 — c Axi). The antidiffusive term can be 
merged with the upwind term by a slight manipula- 
tion to obtain a purely upwind space-marching KNM. 
First, we rewrite equation (89) as 


(®)F(x! +Ax 1) x 2 )> = <*,F> 

-(* r “ lT '-£;( F+i r PcEF )) 

+ ^(*' E2 Sf} +0 ( 4 -'> < 91 > 

With the obvious relations 


<¥,F) = ^f(i + ^P ce ^ 



+ O (Axf) 


and the definition 

fcE = F (l + — 2 ^Ce) ( 92 ) 

equation (91) takes on the simple form 
(*,F(X! + Ax 1 ,x 2 )> = (®,/ce ( X 1 > x 2 - cAii)) + O (Axf) 

(93) 

The velocity distribution fcE is the Chapman- 
Enskog distribution function, and because of the pos- 
itive sign in front of the second term in the brackets 
it is antidiffusive. Thus, the antidiffusion term in 
equation (90) appears in the form of a perturbation 
term in equation (92). The space-marching scheme 
in equation (93) now takes on a purely upwind form. 
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more familiar form 


(p« l )” +1 = J J J "l^CE { x l n ’ x 2 J -cAnj dv 1 

(94) 

(p + P «?)” +1 = j J J v lfCE ( x l n i x, lj - cAzi j dv i <fv 2 dl 

(95) 

(P“l“ 2)" +1 = J j J v l v 2fcE - cAijj dv 1 dv 2 dl 

(96) 


(peu! +P«l)" +1 = ( — ^—PU! + — ■ “ X 
J \ 7 - 1 2 


1,1 /+ T 


x l n , x 2 , - c Ai 


^ dv\ dv 2 < 


The only difference between the first- and second- 
order accurate schemes is that the first-order scheme 
uses the Maxwellian distribution and the second- 
order scheme uses the Chapman-Enskog distribution. 

The integration with respect to / in equations (94) 
to (97) can be performed in closed form. Equa- 
tions (94) to (97) then become 

(P«l )” +1 = / / *>1/CE (zin^ - c Az x ) dvi dv 2 

(98) 

/ OO r OO 

/ V l/CE Pl»,l2y -cAi! j dvjcit^ 

(99) 


(p«lU 2 ) 


+1 = / / «1«2/ce ( I 1„ I 2 > -c Ax x j 

J —oo J —oo ' 

( 100 ) 

/ oo /• oo 

/ 

-oo J — oo 

X '“lJcE ^ x ln > X 2j -^Alij 
X /o ^12^ — c All ^ dttj dv 2 

/ OO /'OO o 

/ T 

-00^—00 

x /CE (*l„.*2y - c Ai x j duj dt; 2 


where 


/CE = E ^1 + ^--Pce 

JcE = f (l+^^CE 


F = Contracted local Maxwellian distribution 
P ..... T {vi ~ urf + (v 2 - u 2 j 2 } 


2t tRT 


exp - 




,~ d P 
+ C dx 2 


)[ 2 ~ 


("1 “ u l) 2 + ( v 2 - u 2 ) 2 


_ du 2 \ v 2 - u 2 
+ c- — — 


( dp 


(t»l - til) 2 + (tt2 - «2) 2 , 

V dx\ 

dx 2 ) 

2 RT 

(K 


(J>L + 

clL) 

j (til ~ ui) 2 + (v 2 - u 2 ) 2 

Vdxi 

dx 2 J 

2 RT 

(^1 4 


t>i - tii / du 2 „du 2 \ v 2 - u 2 

\ dxi 

dx 2 J 

RT Vdx x C dx 2 ) RT 

1 / dp 

. dp 

\ (tij - Ul ) 2 + (v 2 - u 2 ) 2 

P 

dx 2 

) 2 RT 

(103) 


The partial derivatives of p, ui, u 2 , and p with re- 
spect to x\ in equations (102) and (103) are replaced 
in terms of their x 2 derivatives by using the equations 
of appendix D. The x 2 derivative can be evaluated at 
mesh points by using, for example, central differenc- 
ing. The numerical solution obtained with the KNM 
scheme will then be second-order accurate both in x\ 
and x 2 , provided the integration with respect to v\ 
and v 2 is done sufficiently accurately. 

Computed Results 

In this section, we consider the application of 
the kinetic numerical method (KNM) to two test 
cases. The time-marching KNM developed previ- 
ously for the unsteady Euler equations is applied to 
the one-dimensional shock-tube problem. The space- 
marching KNM is then applied to the shock-reflection 
problem. 

Application of KNM to Shock-Tube Problem 

For the application of the KNM to a one- 
dimensional shock-tube problem, we assume that the 
pressure p, temperature T, density p, and velocity u 
at time t = 0 sec are given by 

^HP = PHP = 5.0 


x . 

t lp 

n H P = U LP = 0 
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where HP and LP denote, respectively, the values 
corresponding to the high-pressure and low-pressure 
sides of the shock tube. The initial value of temper- 
ature is assumed to be 300 K and the gas constant 
per unit mass is R — 287 J/K. Results of computa- 
tions are shown in figures 1 to 6. Distributions and 
interpolation schemes assumed in the computations 
are summarized in the following table. 


Fig. 

Distribution 

Interpolation 

scheme 

Courant 

number 

J 

sec 

b 2 

Maxwellian 

Linear 

0.95 

101 

0.00083 

3,4 

Chapman-Enskog 

Eq. (63) 

.95 

101 

.00041 

5, 6 

Chapman-Enskog 

Eqs. (63) 
and (65) 

.95 

101 

.00041 


In this table J is the number of mesh points used and 
t j 4 is the time up to which the solution is advanced. 
The Neumann boundary condition has been used 
for all computations. In the KNM formulation, 
boundary conditions are required to determine the 
values of F ( t n , Xj — v At, v) whenever Xj — v At is 
outside the computational domain xy i < x < xyjj. 
For the Neumann boundary condition at xy — xy i 
and xy = xy jj, the values of F [xj — v At) outside 
the computational domain are given by 


F( Xj 


( F(x 1iL ,v,I ) (xj - v At < xy tL ) 

l F (zi,U,t>,/) (xj - V At > Xy^jj) 


Plots of field variables against x presented in fig- 
ures 1 and 2 show no pre- and post-shock oscilla- 
tions. This is expected because of the TVNI prop- 
erty of the first-order accurate KNM. However, the 
shock and contact surfaces are badly smeared. This 
smearing confirms that the first-order KNM contains 
a large amount of numerical viscosity. The results 
of computations with the Chapman-Enskog distri- 
bution are shown in figures 3 and 4. The shock is 
much sharper, and this confirms the analysis of the 
preceding section. It is shown that the Chapman- 
Enskog distribution reduces the numerical viscosity 
in the KNM considerably and the numerical viscos- 
ity goes like 0(At 2 ). The sharpness of the shock, 
therefore, depends on At. The velocity and tem- 
perature, however, show marked fluctuations in the 
neighborhood of the shock as well as in the proxim- 
ity of the contact surface. Use of the interpolation 
scheme (eqs. (63) and (65)) makes the KNM satisfy 
the TVNI property. This modified second-order ac- 
curate KNM is expected to do much better, and the 
results shown in figures 5 and 6 do indeed confirm 


this. The oscillations near the shock and the contact 
surface are considerably suppressed. However, the 
density plots of figures 3 and 5 reveal that the shock 
has been smeared somewhat as a consequence of the 
TVNI modification compared with the case when the 
interpolation is done with equation (63). 

As mentioned before, the Chapman-Enskog dis- 
tribution is not uniformly valid in u-space and in fact 
becomes negative for v > O (At -1 / 3 ) . Thus, the dis- 
tribution assumed in the computations is not right 
for fast molecules which penetrate through the dis- 
continuity. It is here that work of Narasimha (ref. 11) 
on uniformly asymptotic expansion for / may help a 
great deal. In particular, the precursor distribution 
moving downstream in the neighborhood of the shock 
is likely to improve the KNM even more. 

Application of KNM to Shock-Reflection 

Problem 

The second-order accurate space-marching KNM 
developed previously is applied to the case of an 
oblique shock reflecting from a flat plate. Figure 7 
shows a Cartesian grid system with 61 points along 
the Xi-axis and 21 points along the V2-axis. The 
flat plate is on the Aq-axis and the computational 
domain is 

X 1,L <xy< x \,u 
X 2,L < x 2< X 2,U 

An oblique shock with Mach number My =2.9 and 
deflection angle 8 — 11° enters the computational 
domain at the point (xy t L, x 2 ,[/)• The shock, after 
undergoing reflection at the bottom flat plate, exits 
through the right-hand boundary of the domain. 
The initial conditions for the left-hand boundary are 
obtained from the uniform flow parallel to the Xy- 
axis at all mesh points except (aq £, £2,t/)- At this 
point the flow corresponds to conditions downstream 
of the oblique shock. 

To advance the solution from line n to n + 1, 
a scheme for numerical quadrature and a bound- 
ary strategy are required. The quadrature in equa- 
tions (98) to (101) is performed with the product rule 

h{x 1 ,x 2 )dxidx2 = YT”^h (zij.z 2> ) 

l j 

where w. L and Wj are the weights for the one- 
dimensional Simpson quadrature scheme. While per- 
forming the numerical quadrature, we need to evalu- 
ate the Chapman-Enskog distribution function / q e 
at various mesh points. For some values of c, the 
point X 2 - — eAaq may be out of the domain. Con- 
sequently, boundary conditions are then required to 
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calculate /q e \xi n ,X2j — cAxij. Two possibilities 

arise: 22 — c At < X2 or 22 — c At > 

For the shock- reflection problem, the flow tan- 
gency condition has to be imposed at the lower 
boundary (which is a solid wall), and the Neumann 
condition can be used at the upper boundary (which 
is fictitious). The details are as follows. Let x ' 2 
be the position from which fluid particles arrive at 
X2 during the time-like interval Ax\. Obviously, 
x ' 2 — X2 - cAx\ if 22, L < 22 cAx\ < X2 jj. If, 

however, 22 - cAxi > 22,(7, then set x 2 = X 2 ,U ■ 
This ensures that the Neumann boundary condition 
is imposed on the upper boundary. 

For positive values of c, the point 22 — c Ax^ can 
be below the lower wall, that is, 22 — cAxi can be 
less than 22 l- In such a case, a specular reflection 
model is used to satisfy the flow tangency condition. 
The value of x ' 2 can be obtained with 


Finally, to evaluate /qe \ x i n -, x 2- — c AxiJ when- 
ever X2- — c Ax\ is not coincident with a mesh point, 
an interpolation scheme is essential. The flux limit- 
ing method described previously can be used to make 
the KNM a TVD method. 

The results of the computation are shown in 
figures 8 to 12. The pressure contours obtained 
with the flux limiting interpolation are depicted in 
figure 8. Figures 9 to 12 show the pressure plots 
at two stations. Also shown on these plots are the 
results with second- and third-order TVD methods 
of Chakravarthy and Osher (ref. 12). The position 
of the reflected shock given by the KNM agrees 
better with the position given by the third-order 
TVD method than with that given by the second- 
order TVD method. 

Review and Discussion 


x 2 = 22 L + c (A21 — Axi) 


Axj = 


x 2 ~ X 2,L 
c 


The particle velocity after the reflection at the wall 
is — c. Physically speaking, the gas-surface model 
consists of the fluid particle starting at x' 2 , traveling 
for A21 — A21 at — c, and then traveling up to 22 for 
the interval A21 with reversed velocity. For each pair 
of values 22 and c, we determine x ' 2 and c' and then 
use these values in calculating /qe ( 2i n , 2^, v\,c') 
to be further used in the numerical integration. 

The numerical integration in equations (98) to 
(101) generates the conservative variables 


U 4 - (pe + p)ui = 


U 1 = pui 

U 2 =p + pu\ 

f/ 3 = puiu 2 
/ 

7 


;P 


pu 


U 1 


7- V 2 
The field variables p, iq, «2, and P are given by 

. 11/2 


7 - 1 


+ 


Ui = 


( 7 A) -^1(2 c,t/ 4 -c|) 




7 
Ui 

p = — 
u 1 

_ ^3 

U 2 Ui 

p — U2 - UlUi 


A new upwind method, called the kinetic numeri- 
cal method (KNM), for obtaining the numerical solu- 
tion of the Euler equations has been presented. The 
method exploits the well-known fact that the mo- 
ments of the Boltzmann equation of the kinetic the- 
ory of gases are the Euler equations when the distri- 
bution function is Maxwellian. Thus, every method 
for the solution of the Boltzmann equation can be 
mapped to a method for the solution of the Eu- 
ler equations provided the distribution function is 
Maxwellian. The mapped-method strategy has been 
used to develop the KNM, which is explicit, is up- 
winding, satisfies the entropy condition, and is highly 
vectorizable. Further, it can be made total variation 
diminishing (TVD) by using a suitable interpolation 
strategy for evaluating the Maxwellian distribution. 
A unique aspect of the present work is the use of the 
antidiffusive Chapman-Enskog distribution in devel- 
oping a second-order accurate KNM. An immediate 
consequence of the use of the Chapman-Enskog dis- 
tribution is that by modifying slightly the expressions 
for r and q to include physical viscous stress and heat 
flux, the KNM can be applied to obtain the solution 
of the Navier-Stokes equations. 

An important difference between the KNM and 
other methods is the treatment of the boundary con- 
ditions. On any surface (bounding the computational 
domain) where boundary conditions are imposed, the 
velocity distribution function / can be regarded as 
consisting of / + and standing respectively for the 
particles moving towards and away from the surface. 
The particles corresponding to / + carry information 
out of the domain, whereas those corresponding to 
f~ cause an inflow of information from outside. It is 
then possible to treat the boundary conditions in a 
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physically meaningful way. For example, in the ap- 
plication of the KNM to the two-dimensional shock- 
reflection problem, particles reflect with reversed nor- 
mal velocity after colliding with the solid wall. The 
distribution f + is then the velocity distribution due 
to incident particles, whereas / - corresponds to the 
reflected particles. The flow tangency condition is 
then satisfied by simply requiring that 

/“ (v n ,v t ) = / + (-V n ,V t ) On < 0 ) 

Another unique feature of the KNM is that the 
solution at every field point X2 is due to the particles 
directly arriving from any other point x' 2 - The 
solution is also due to particles starting from x' 2 
undergoing interaction with the solid surface and 
then arriving at X2 with different particle velocities. 


In terms of / + and f~, we can state that the solution 
at X2 is influenced by both / + and /“. For field 
points close to the solid boundary, the solution will 
be strongly influenced by particles reflected from 
the boundary. These two kinds of contributions to 
the solution are physically meaningful and allow' the 
KNM to treat the boundary conditions in a way 
that is not found in many other finite-difference 
schemes. In view of the above advantages of KNM, it 
is believed that the approach of developing mapped 
methods will lead to the construction of numerical 
schemes for the Euler and Navier-Stokes equations 
not known heretofore. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
July 10, 1986 
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Appendix A 

The Connection Between the Boltzmann Equation and the Euler Equations 

The Boltzmann equation for / is 


0 / , df 
aT + v 'fc 


= AU) 


The collision term J has collisional invariants 1, v, and /+ tj-, that is, 


(Al) 


ff ^(v,I)J(fJ)DvdI = 0 (A2) 

2 

where -0 = 1, v, and I + \ and Dv = dv\ dv 2 dv 3. Equation (A2) is a direct consequence of the 
conservation of mass, momentum, and energy during a collision. Substituting equation (Al) 
into (A2) yields 

fh {v ' i) (% + ' , -%) Dvdi=o (a3> 

The field variables density p = p(£,x), fluid velocity u = u(f,x), and specific internal energy 
e = e(t,x) are related to / through the following moment equations: 

P = II f{t,x,v)DvdI | 

pu = 

pe - 

The viscous stress tensor r and the heat flux vector q are defined by 


//v/(*,X,v) 

if 


Dv dl 


(A4) 


/+ y I /(f,x,v) Dv dl 


P&ij T ij 

= 11 C{Cjf Dv dl 

(A5) 

*"Jf 

(l+ C -pJ Ci fDvdI 

(A6) 


where q is the peculiar velocity v % — Uj, b,- is the Kronecker delta, and p is the static pressure. 
Using equations (A4), (A5), and (A6) we have 


r 

1 

df ^ , T 

d 

P 

1 

V 

r , v 2 

-bf-Dvdl = 
dt 

dt 

pu 


L' + T- 


,P e . 


and 



/ + S 


yj 


in- — — Dv dl 

3 ° X J 



v j 



f Dv dl 
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Now, 


J j v^vj f D\ dl = JJ (c{Cj + c^Uj + CjU t + f Dv dl 

P^ij ~ij “I" P u i u j 


and 


*//K 


Ujf Dv dl 


Cjf Dv dl 


—peuj + JJ + 

=peuj - qj + (p6 { j - T^ 'j 


(v - u) 2 u 2 

2 + T + c%Ui 


Cj Dv dl 


u; 


Equations (A4) therefore yield 


dp d 
dt dx 


d_ 

dt 


: ( P u i ) = 

i 

(p u i) + JJr, + Pkj) = 


ij 


d „ 


d_ 

dt 


'j "xj 

d ( \ d 


(P e ) + ~ ( P eu j + P u j) = faT. (<?i + u i T ij) 


(A7) 


The above moment equations reduce to the Euler equations if 

r t j = 0 and qj = 0 

In general, for an arbitrary /, the moments Tjy and qj are not zero. However, if / is the 
Maxwellian distribution defined by 


f = F = 


P 

(2? tRT) 3 / 2 


exp 


(y ~ u ) 2 exp (—/// 0 ) 
2 RT 1 0 


(A 8 ) 


then Tjj and qj vanish. Here R is the gas constant per unit mass, T is the temperature, and I 0 
is the internal energy due to nontranslational degrees of freedom. The value I 0 is related to T 
through the relation 


Io = 


( 2 + Df) - 1 D f 
2(7-1) 


RT 


where Dj- is the number of translational degrees of freedom and 7 is the ratio of the specific 
heats. When / is a Maxwellian distribution, 


JJc^fDv d, = jf c 2 c zf Dv dl = JJ c%c\fDvdI = 0 
JJ c\fDvdI — JJ c\fDvdI = J J c|/ D\ dl = pRT 
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The normal stress tensor components rn, T22, and T33 can be absorbed in the definition of p , 
and equation (A5) then yields 


V = 



c\f D\ dl = pRT 


(A9) 


22 


Appendix B 

Elimination of Time Derivatives 

The one-dimensional Euler equations are 

d 3 

at w + ai ( '"‘ ) = 0 

d d 

rt (pu) + -^(p + ru 2 ) = o 
d . , d , 

ar pe ) + fa(f )ue + P u ) = 0 

where 

RT u 2 p v? 

e = 1 = — - — — — - H 

7-1 2 p(7 - 1) 2 

The pressure p is related to e through 


(Bl) 

(B2) 

(B3) 


(B4) 


P = p{l~ 1) e - 


The objective is to replace 


§i M ’ f<' ra) ' 1^)' |lW’ and 

in terms of the space derivatives. Taking p first, we have 


3 2 dd. d d . 3 2 . 2 \ 

a?*"* = -atai ( '" ,) = -&a+" ) = a? (p + ' m > 


Taking pu next, we have 


3 . . 3 2 \ 

d 2 dd 2 \ 

a? (p “ ) = -toa (, ’ + p “ ) 


Now, the energy equation can be written as 


3 / 7 — 1 2 \ 5 / 7 + 1 3 

^(p+— pu)+^bpu + — pu 1=0 


(B5) 


(B6) 


(B7) 

(B8) 


(B9) 


which is obtained by eliminating e in terms of p in equation (B3) through the use of equation (B4). 
Hence, 


d 2\ & 

m {p+im} = dt 


P + 


7-1 


-pu 


d_ 

dx 


^7 pu + 


1 ~ 1 

2 


3 — 73 
2 dt 


+ 


pu° + 


(pu 2 ) 


3-7 d , 2 

—3+“ 


(BIO) 
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as a result of equation (B9). For the time derivative of pu 2 we have 


d_ 

dt 


{pu 2 ) = u 2 ^(p) + 2 Pu^~ t {u) 


dt 

_ d o 3 , . 

= 2 u-(,m) - u -(„) 

d 2 2 ^ t \ 

= _ 2 „_ _ (H 


where equations (Bl) and (B2) have been used. Equation (BIO) then yields 


u l 1 ~ 1 3 — 1 2 & t \ 

' “ *"“• ' “ ^~(H 


^( P + P » 2 ) = -^I 1P « + 


3x 


-( 3 -7)w^(P + P« 2 ) 


Slight manipulation gives 


d 2 \ 3 

Tt {p + pu ) = "ai 


'-v — 1 o o 

7pu H — pu + (3 - -y)pu + (3 - -y)pu 


3-7 3 

— 


+ (3 - 7) (p + pu 2 - pu 2 ) — (u) 


(9 (9 

= - ^ (3P« + + (3 - q)p— (u) 

Note that 

JJv i FdvdI = JJ j(w — u) 3 + 3(u - u) 2 u + 3(v - u)u 2 4 - u 3 ) Fdvdl 

o 

= 3 pu + pu 

Equation (Bll) can then be written as 

d , 9, d f f o , „ ,3 


-(p + pu 2 ) = --JJ v * FdvdI + {3- 7 )p- ( 


Combining equations (B12) and (B8) gives 


3 2 . . d 2 o d 

S 2(H = ^(3p« + P » )-^ 

f v 3 Fdvdl - 

dx 


d 2 [ 

dx 2 J 


(3-l)p^W 


Finally, let us consider pe. Its second-order time derivative is given by 

3 2 3 3 

For the time derivative of pue + pu we have 


3 3 I 

— (f,eu + pu)= — lu 


Pe + p(7~ 1) ( e- — 
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(Bll) 


(B12) 


(B13) 


(B14) 


where p has been eliminated in terms of e with equation (B5). After slight simplification we obtain 

3 , 3 7-13.0 

-( pen + = -{pent ) ) 

3 / pu\ 7 — 13 p 3 u 3 


= d ( oe PJi\ ^~ ld P u 

dt \ p ) 2 dt p 2 

[ 3 3 3 


„3b_j) u 2 a w + h _ i)u 3a w 

which is a relation containing first-order time derivatives of p, pu, and pe. Grouping terms together 
we have 


3 3 f 3(7 — 1 ) 2 1 3 

— (peu + pu) = 7 U— (pe) + 7 e u — (pu) 

o 1 3 

+ [(7-l)u -7uej— (p) 

Using equations (Bl), (B2), and (B3) to eliminate time derivatives, we obtain 

3 3 [ 3(7 - 1 ) 2 1 3 2 \ 

— (peu + pu) = - 7 U— (peu + pu) - 7 e r — (p + pu^) 

o 1 3 

- (7 — l)u -quej — (pu) 

Substitution for e in terms of p through the use of equation (B4) yields 


(B15) 


3 3/7 pu 3 

+ = p»+ — 


3 2 x 

x^(p + pu)-[ 


7 P . 3 - 27 2 

7— lp + 2 
7P W 13,, 


7-2 o 7 pu 3 , , 


We can replace each term on the right-hand side of equation (B16) using the identity 


4< b >=I (ab) - b I< a > 


Equation (B16) then reduces to 


(B16) 


3, , 3/ 7 2 274 7 p 2 7 2 3 — 27 2 

— (peu + pu) = - — -pu 2 + ^pu 4 H 7 — + 7 PU + — 7—^ Pu 

dt 3x \ 7 — 1 2 7 — 1 p 7 -I 2 

3 — 27 4 7 — 2 4 7 2 \ 

+ -^pu 4 + — /m 4 - — pu 2 j 




J-!t ,, d It- 2 3 
+ ^-» ) + ^(— “ 


7 pu 

P ( 7 - 1) 


3 57 — 3 2 7 P 2 1 4 /7 2 7 3 \ 3u 

aT + UT7 + +" + (uT f ’" + 2' m ) 


+ (p + ' , ” 2) £(uu + 


3-27 2 \ , 3 [7-2 3 7PU 

— " “ “TUT) 
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Finally, grouping terms containing j^u and ^ we obtain 


d , , d 

— {peu + pu) =- — 


Note that 


dt 


dx 


5^-3 2 


2 ( 7 -l) 


pu + 


IP, 1 4 

r— + -pu 1 

7 ~ 1 P 2 


+ £ w 


+ 7jPU 3 + (p + pu 2 ){ 3 - 27 )u + ^(7 - l)pu 3 
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P7 , 
7 - 1 


+ 


/ . 2\ 1 2 7 

(p + ptT)- 7 - pu 


d_ 

dx 


57-3 2 


+ 


2(7-1) 
IP d (p 


pu + 


7 - 1 

7 P 2 , 1 4 

7-1 P 2 PU 


A (V 

7 - lj \P, 


5 


+ (3-7)pu— (u) 


7 — 1 dx \p 



1 + 


v 2 F dv dl = JJ Iv 2 FdvdI + ^JJ v 4 FdvdI 

— Io{p + pu 2 ) + ^ JJ [(u - u ) 4 + 4(v - u) 3 u + 6(u - it) 2 


+ 4(w — u)u 3 + u 4 


Fdv dl 


-=7} ~ P ( p + pu ' 0 + 5 (t - + 6pu2 + pu4 ) 


3-7 p 

2(7 


57-3 _ 2 


2(7-1) 

Hence, equation (B17) can be written as 


d d 

-(peu + pu)^-- 


pu“ + 


7 P , 1 4 

q + -^pu 

7-1 P 2 



/ + — v 2 F dvdl 
“ I 


+ (3 - 7 )pu—u 
dx 


+ 


IP d (p 


7-1 dx \p, 

The second derivative of p with respect to t is thus given by 


dt * (pe) ~ dx 2 



dx [7 — 1 dx 


\ 2 1 

d 

du 

v 2 F dv dl 

dx 




(B17) 


2 


(B18) 


(B19) 
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Appendix C 

Chapman-Enskog Analysis 

Consider the model Boltzmann equation 


3 / , !>i 


where 


F = J-(~ ) exp \-/3(v - u) 2 - I/I a 

1 0 \7T/ t 


7 . 3-, 3 — "y 1 ^ 

° 2( 7 -l) 4(7 — 1) /? J 

and a is a time constant characteristic of collisions. Assume that / = where $ is a function to 
be determined. Substitution in equation (Cl) yields 


_ , $ 1 (dF dF\ F (d* d$ 

®~ 1 ~aF\lH +V fa)-a{lk +V lte 


Formal solution of equation (C4) can be written as 


_ 1 — (F/a)[(d$/dt) + v(d$/dx)] 
~ 1 + (1 faF)[{dF/dt) + v{dF/dx)} 


F (d$ d§\ $ 1 (dF dF 

a\~dt +V l)x ) << aF\l)t +V lte 


then approximate solution for $ is given by 

1 + (1 /aF)[(dF/dt) + v(dF/dx)\ 

Further, if ^ ) << 1, then $ « 1 - + vff ) 

The Chapman-Enskog solution is 


, ^ r 1 (dF OF 

JCE aF\dt dx 


The next step in the Chapman-Enskog analysis consists of replacing ^ with space derivatives 
determined with the Euler equations. For this purpose note that 


1 (dF dF\ ( d n 

F(ar + "fe) = U + t W lnP 


In F = In p + iln In 7r — f3(v — u ) 2 -\nl 0 — ^- 
Z Z 1q 
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Hence 


1 (dF dF 

Pce = f (^F + V a^ 


_ 1 (dp dp 
- ~ P \Tt +v d~x 


2 (3[v — u ) 


1 

20 “ 
du du 

m +v di 


+ t 72 

rt 


d(3_ 

dt 

I 

T 0 


d/3 

v !Tc 

dlo 

dt 


dlo 

dx 


(C7) 


From equation (C3) it follows that 


Io\dt dx ) I3\dt dx ) 


Equation (C7) therefore becomes 




P 

+ 


, dt 

A _ 1 

.2/3 B 


I + 20(v 
( v — u ) 2 


( du du 
t m +v d~x 


9/3 

Tt 



(C8) 


where 


Now, from the equation of continuity, 


B = 


3-7 
4(7 - 1 ) 


dp dp du 

dt U dx P dx 


From the momentum equation, 

du 1 dp du 

dt p dx dx 

_ 1 d ( p \ du 

p dx \2(3 ) dx 
1 dp 19/3 du 
2p(3 dx 2/3 2 dx U dx 


Thus, 



du du du 1 9/3 1 dp 

dt JrV dx V U ' dx + 2/3 2 dx 2p(3 dx , 


(C9) 


Substitution of equations (C9) into (C8) yields 


P C e = 2/3 {v - uf 


du v — u d/3 f d(3 9/3 

dx /3 dx ^ \ dt dx 


3 7 2 

— — (v — u) 

1.2/3 B K ’ 


(CIO) 


Now 


If 


remains to be eliminated. The energy equation gives 


9e 9e 9 . 

p ai +l,u Tx + Vx (pu) 
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Substituting for p and e in terms of P by using equation (B4) and the equation of state, we obtain 


1 d(3 du 

2(7 - l)/3 2 dt + U dt 
p du dp 


+ pu 


1 d(3 du 

+ u 


2(7 — 1 )/3 2 dx dx 


Using the momentum equation to remove we obtain 


or 


PP d P , , . R du 


d0 d/3 dp du 


Substituting «f + »M from equation (Cll) into (CIO) gives 


P CE = mv ~ u) 2 - 1 


du v — u dp 
dx + P dx 


3 I , ,2 

12 . 




3 7~5 , /, ^ af .. .02 (7-1)^ 


+ 


+ (3 - 7 )P(v - up 
5(u — u) I(v — u) 


B 


du 

dx 


( v — uf 


2 p B 

Substituting for B in equation (C12) gives the following result 


dp 

dx 


FCV = | —^7, “ + (3 - 7 )P(v ~ uf - A \~tpl 

+ I T '-I{V -U)-{V~ uf 


b(v — u) 4(7 - 1) 


du 

dx 

dP 

dx 


2 P 3-7 

Note that Pqe * s a polynomial in v — u and I. The Chapman-Enskog distribution 

PCE 


Ice = F 


a 


(Cll) 


(C12) 


(C13) 


(C14) 


is thus a perturbation over the local Maxwellian distribution F. Some important properties of /ce 
are worth noting. First, observe that 

JJ F dvdl — p, JJvFdvdI = pu, j J v 2 F dv dl = p + pu 2 , 

JJ (v — ufFdvdl = p, J J (v — u)F dv dl — 0, JJ IFdvdI = pI 0 , 

JJ I 2 F dvdl — 2pl 2 , JJ [v - ufFdvdl = 0, JJ (v - u) 4 F dvdl — 

JJ {v - ufFdvdl = 


29 


With these definite integrals the following results are obtained: 



// 


1 

v 


I + \ 


(/ce “ F)dvdl = 0 
ff ( v ~ 0 2 /ce dvdl = p- r CE 


/ + 


(w — a ) 2 


{v-u) fcE dv dl = q CE = 


7 V d (P 


7—1 a dx \p 


where 


TCE 


3 — 7 du 

PR~ 

a ox 


Observing that 

dx \p) 2/3 2 dx 

We can write the Chapman-Enskog distribution terms of r£ E and qc E as 


fCE = F 


T CE p _ QCE p 
P T pV2RT q . 


where the polynomials P T and P q are given by 


where 


Pr(v,I) 


3^-5 r 2 _ 4 ( 7 ~ 1) 2 I 
2(3 - 7 ) (3 - 7) 2 2RT 


P q {v,I) 


2(7-1) 

7 



4 (7 ~f) r I 
3-7 2RT 


(2i?T) 1 / 2 


The kinetic numerical method involves integrals of fc E with respect to v and I. The 
momentum transport involves 


Jf fcEdvdl and // v fCE 

whereas the energy transport involves 

JJ I f CE dvdl and H^fcE 

The dependence of I can be integrated. For 1-, v-, and u 2 -transports, we then have 
fcE = I fcE dl 


dv dl 


dv dl 


= F 


T CE 


( c2 -0- 


QCE 2(7 - 1) ( c3 _ 3^ 


p(2i?T) 1 / 2 7 


where 


F = 


tfnRT ) 1 / 2 


exp (— C 2 ) 


(CIS) 

(C16) 

(C17) 

(CIS) 

(C19) 

(C20) 

mass and 


(C21) 
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For /-transport 


fCE = j IfCE dI 

_ 2(2-1 ) 9CE (c 3 ~ V)} (C22) 

1 p{2RT) 1 / 2 V 2 /] 

not depend on / but do depend on v. The second- 


[:»] 

= J [J] 7ce (x-vAt)dv 


(C23) 

(pe ) n+1 = J 

V 2 - 

Jce{x — vM) + y fCE(x - v At) 

dv 

(C24) 

To obtain <jce an d in equations 

(C23) and (C24), a is equal to 




= F<1- 


T CE 


c 2 


1 ' 


2(3 - 7 ) J 


The contracted distributions fcE> f CE> and F do 
order accurate KNM (eq. (58)) then reduces to 



Appendix D 
Euler Equations 

The two-dimensional Euler equations for a steady flow can be written in the following form: 


dp du\ dp du 2 

Ul dxi + ^dxi U2 dx2 P dx2 

(Dl) 

dui 1 dp du2 

Ul dx 1 + p dx 1 U2 dx 2 

(D2) 

dv,2 1 dp du2 

3xi pdx 2 dx 2 

(D3) 

dp odui dp 2 dv,2 

ui--— + pa - — - -u 2 ^ pa - — 

3xi axj dx 2 3x2 

(D4) 


where 

a 2 = i RT 

The only difference between the above equations and the usual form is that all the space derivatives 
with respect to x\ are on the left-hand side and all the space derivatives with respect to X 2 are on 
the right-hand side. 

Solving equations (D2) and (D4) for f^rj- and we obtain 


du i U\U 2 dui a 2 du 2 U2 dp 

dx\ u 2 — a 2 dx 2 u 2 — a 2 <9x2 p (u 2 — a 2 ) <9x2 

dp _ pa 2 U2 dui pa 2 u\ du2 U\U2 dp 
dx i u 2 — a 2 dx 2 u 2 — a 2 dx 2 u 2 — a 2 dx 2 

Substituting for and * n equation (Dl) and simplifying, we get 


(D5) 

(D6) 


dp 

dxi 


U 2 dp pU 2 du\ 
«i dx 2 uj - a 2 <9x2 


pu\ du2 
u 2 — a 2 dx 2 


U 2 dp_ 

ui ( u\ - a 2 ) dx 2 


(D7) 


The results obtained can be written as the matrix equation 


du\ 
dx 1 


dU2 

dxi 

dp 

. dxi . 



ui 

0 

0 

0 


PU 2 

pu 1 


U2 -1 

u{ - a 2 

uj — a 2 

Ul 

(u 2 — a 2 ) 

“1«2 

a 2 


a 2 

Uj — a 2 

u 2 — a 2 

p{u\-a 2 ) 

0 



1 

ui 


pu 1 

pa 2 u 2 

pa 2 u 1 


ui U 2 

uj - a 2 

2 9 

— a z 


u 2 — a 2 


If' 

ox 2 
du\ 

dx2 

du2 

dx 2 
dp 

L dX2 . 


(D8) 


With equation (D8), the space derivatives with respect to x\ can be expressed in terms of the 
space derivatives with respect to X 2 - The Chapman-Enskog polynomials Pce ar, d P CE defined by 
equations (102) and (103) can thus be determined only in terms of v±, V 2 , p, «i, 112 , p, and the 
space derivatives with respect to X 2 . These space derivatives can be determined by a second-order 
accurate differencing scheme such as a central difference formula. 
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for shock-tube problem with Maxwellian ansatz. 


A P 



shock-tube problem with Maxwellian an: 



1 .2 .3 .4 .3 . 6.7 .8 .9 1.0 


Figure 5. Computed and exact p and u profiles for shock-tube problem with Chapman-Enskog ansatz with 
TVD modification. J = 101; £4 = 0.00041 sec. 





Figure 6. Computed and exact p and T profiles for shock-tube problem with Chapman-Enskog ansatz with 
TVD modification. J = 101; 1 4 = 0.00041 sec. 






Figure 7. Computational grid with 61 x 21 points for reflection of oblique shock by a flat plate. 



Contour from .9000 to 4.1000 
Contour interval is .10000 

Figure 8. Pressure contours for shock-reflection problem. 
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2 


o Second-order 
TVD method 



Figure 9. Pressure profiles computed with KNM and second-order TVD method for shock-reflection problem 
with 61 X 21 grid at x% = 0.0. 





Figure 11. Pressure profiles computed with KNM and third-order TVD method for shock- reflection problem 
with 61 x 21 grid at X 2 = 0.0. 



Figure 12. Pressure profiles computed with KNM and third-order TVD method for shock-reflection problem 
with 61 x 21 grid at X 2 = 0.5. 
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